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SEQUENTIAL  ARRAY  ALGEBRA  USING  DIRECT  SOLUTION  OF  EIGENVECTORS 


PROBLEM  OF  SEQUENTIAL  ARRAY  EQUATIONS 

The  new  computationally  powerful  array  algebra  technology  unifying 
the  sciences  of  numerical  analysis,  mathematical  statistics  and  modern 
signal  processing  would  become  more  flexible  if  the  problem  of  sequential 
array  observation  equations  could  be  efficiently  solved,  Rauhala  (1974 


p  113,  1976  p  79  ,  1977,  1978,  1979,  1980a,  1980b),  Jancaitis  and  Magee 

(1977),  Snay  (1978).  In  the  illustrative  case  of  three  dimensions  the 


sequential  observation  equations  read 
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where  the  array  multiplications  /£“  —  (Z  -  V) 

'*V*r  /?,  /I4  /7j 


are  defined  as 
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The  last  set  of  observation  equations  consists  of  dot  multiplications, 
i.e.,  discrete  direct  observations  of  parameters  £  so  that  in 
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the  conventional  monolinear  notations  where  X  t  /V-  ^  ^  is 

***!  * 

treated  as  a  long  column  vector  the  design  matrix  would  be  diagonal. 

The  above  observation  equations  result  in  the  normal  equations 


/i  C*i o 

*,  z.  +  -  /?. Z3r  *  z  =  (£,  ,  <» 

^  'VV’j  3 

where  the  dot  multiplications  ^  //  a*e  ^enotec*  ^  • 

We  now  assume  that  the  symmetric  square  matrices  £  are  brought  to 

satisfy  the  following  spectral  decompositions,  for  example  by  using  the 


parameter 

transformations 

of  Buchanan  and  Thomas 

(1968) , 
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Thus  ^  is  the  common  orthonormal  eigenmatrix  of  all  matrices  /JL and  Sj  7" 
are  its  counterparts  of  matrices  Cy  $  /:  /,  •  ••  jb  .  The 

diagonal  matrices  /^t'y  & *  contain  t^ie  eigenvalues  of  matrices  X^d^Cy, 

The  present  paper  is  focused  on  the  computational  solution  of 
equation  (3)  under  the  spectral  decomposition  of  (4).  The  derivational 
part  of  the  solution  is  rather  straight-forward,  i.e.,  premultiplications 

r 

with  71  ,  post  multiplications  with  S  and  the  "back"  multiplications 
with  T^esult  in  the  solution  of  the  diagonal  system  by 

tt  y~r 

X  Z Sr  =  W*  X  u/Jr 
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Now  the  inverse  transformations  with  St  7" result  in  the  solution 
familiar  from  the  filtering  theory  of  signal  processing 


X  =  Xr{//*  X  WSrJ  S  . 


(6) 


In  terms  of  signal  processing can  be  called  "transfer  function". 

In  terms  of  the  general  theory  of  linear  estimators  and  matrix  inverses, 

A 

Rauhala  (1980b),  estimator  X  unbiased  if  all  />.'  •  •  C 

Wj 

For  biased  or  nearly  biased  parameters  .  A  -  •  •  — ^  oo  »  the 

7  4/i/j 

A 

bias,  variances  and  the  norm  of  can  be  minimized  through  the 

pseudo-inverse  solution  simply  by  putting  /tt  :  •  -  O  for 

/»vVj  Wj 

All  of  these  solutions  of  normal  equations  satisfy  the  least  squares 
criteria 


HV,//  v-  //!$//+  •••  //  //  //V*//  -  /f7/S7 . 
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In  several  applications  of  array  algebra  the  dimensions 

A 

of  the  array  ^can  range  several  hundreds  so  that  the  array  solution 
of  millions  of  parameters  is  split  into  the  problems  of  solving  three 
small  orthonormal  eigenmatrices  7*.  After  these  matrices  are 

known  the  array  multiplications  of  equation  (6)  can  be  performed  along 
the  lines  of  the  computer  program  presented  in  (Rauhala,  1980a).  The 
remainder  question  of  this  paper  handles  the  computational  problem  of 
solving  for  matrices  Jj  ~7“. 
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DIRECT  SOLUTION  OF  EIGENVECTORS 

The  computation  of  eigenvalues  Ajof  matrix  A  and  the  corres- 

ponding  eigenvectors  is  presently  dominated  by  iterative  methods  putting 

severe  restrictions  on  the  dimensions  and  conditioning  of  the  matrix. 

Further  the  iterative  solutions  do  not  guarantee  the  orthonormality  of 

matrices  in  -  72t  A  . 

Aft  An  A  A  Art 

In  the  new  direct  approach  of  f  inding  JP^we  split  the  eigen¬ 
value  problem  in  two  separate  parts,  i.e.,  we  assume  that  the  eigenvalue 
i  •  is  known  or  computed  a  priori.  We  are  seeking  direct  solutions  for 

the  corresponding  eigenvectors  %  yP  as  the  non-homogeneous 

A,i  t,n 

solution  of  the  consistent  systems 


X,-  r  O  C8a) 

y,r  *,•  *  °  >  (8w 

where  . 

=  /?  -  Af  Z  .  (9) 

The  solutions  are  found  using  the  general  theory  of  matrix  inverses, 
Rauhala  (1980b),  by 

xt  =  (Z  -  (10a) 

SjT  *  Uf(Z-  A?,/?*)  , 

(10b) 


Vectors  4£*are  arbitrary  and  the  g-inverse  /?t‘  needs  to  satisfy 
the  condition  #4- /9? rff  -  /9j  in  order  to  have  (10a),  (10b)  as  the 


solutions  of  (8a),  (8b). 


Because  by  the  definition  efe? /&*/  ~  / &  •AfX/  =  O 

the  maximum  rank  of  matrix  /?.•  is  fz/*-/  .  We  perform  the  rank 

Aft 

factorization  of  /?4*  as 

(U) 

where  the  submatrix^ has  to  satisfy  the  condition 

This  condition  can  be  derived  by  eliminating  the  ’’independent’’  parameters 

Zf  from  the  system 

rti 

Z,  +  #p  r 

^  2^  t  2^  “ 

by 

z,  *  (*,-*,  z<> 

Substitution  into  the  linearly  dependent  part  of  (12b),  yields 
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The  computational  rule  of  finding  of  (11)  is  simply 
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because  it  satisfies  the  condition  /9f  yielding 
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The  unnormalized  eigenvector  solutions  become  now  from  (10a),  (10b) 
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where  the  last  terms  of  Ut  If,  axe  chosen  to  be  ones.  The  normalized 

Aii  t,n 

eigenvectors  become 
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By  repeating  these  solutions  for  all  eigenvalues  /t  ...  /?  the 
sought  orthonormal i zed  matrices  ^  of  A  JZj.  become 


/>/» 
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[ 


(19) 


For  a  symmetric  matrix  a  *  7Z, , 


OPERATION  COUNT  OF  THE  DIRECT  SOLUTION 

Computation  of  for  a  symmetric  non-sparse  matrix  requires  in 

the  order  of  ft  **  operations  (scalar  multiplications  and  additions) 

for  each  MV)/  s  Z  >  or  totally  /?  operations.  In  our  array  algebra 
'z' 

solution  this  operation  count  is  by  no  means  prohibitive  as  we  are  solving 
for  /V*  s*ifijfl.  ~  parameters.  In  fact,  the  three-dimensional 

array  multiplications  using  the  general  non-sparse  matrices  ^,^7“  in 
equation  (6)  require  the  same  magnitude  of  operations,  Rauhala  (1976, 
1979,  1980a),  Blaha  (1977). 

If  the  spectral  decomposition  of  the  leading  partition 

■»  r  *»•  «• 

of  /?  were  known  as  Z?  A  A.  and  the  same  partition  could  be 
used  for  all(J4k^*then  we  could  perform  the  one-time  multiplication 

^  j .  and  each  fcj/  *  *  would  require 

/  ♦S'*-/**,/  r 

yy-^operations  or  totally  ^  would  require  operations.  This  is  the 

/*/» 

same  magnitude  of  operations  required  for  a  two-dimensional  array  multi¬ 
plication  of  the  type  of  equation  (6). 

In  several  practical  applications  matrix  is  banded  and  only  a 
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few  last  terms  of  are  non-zeroes  so  that  each  ^J^.requires  6 

operations  or  totally  j *  operations  for  7Z.  ,  where  &  is  the  bandwidth 

/?/» 

(usually  6<-S  for  symmetric  matrices).  In  some  practical  experiments 
the  author  performed  the  double  precision  orthonormalization  of  a 
300  X  300  tridiagonal  matrix  in  a  CPU  time  of  a  few  seconds  using  a 
minicomputer. 

APPLICATIONS 

The  above  array  solutions  were  used  for  simulations  of  non- 
separable  filters  of  finite  element  solution  of  regularly  gridded  data. 

Using  these  filters  or  impulse  responses  a  rigorous  least  squares  solution 
of  601  X  1201  >  720  000  nodes  was  convolved  in  a  CPU  time  of  less  than 
one  minute  and  using  less  than  30  K  bytes  of  the  minicomputer  core  space. 

For  the  non-stat ionary  case  of  irregular  gridded  data  the  above 

derived  array  solution  (6)  removes  some  restrictions  of  the  one-batch 

array  equations.  For  example  in  digital  terrain,  geoid,  gravity  anomaly 

etc.  modeling  using  the  method  of  array  algebra  finite  elements  the 

observed  nodes  are  allowed  to  have  completely  arbitrary  locations  and 

a  priori  weights.  Simultaneously  the  operators  R,fSt’7m  can  be  brought 

to  exhibit  a  structure  of  generalized  fast  transforms,  (Rauhala  (1980a), 

so  that  are  never  explicitely  computed  (requiring  no  core  space) 

and  multicplication  ft  hf  requires  less  than  operations,  i.e. , 

/>/>  n,  t 

the  total  solution  of  /V  r  -  /J3  parameters  requires  the  magnitude  of 

A/ operations. 

The  above  very  fast  array  solution  can  exhibit  such  general 
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properties  which  seem  to  be  at  odds  with  the  restrictive  nature  of  each 
sequential  array  batch:  For  example,  the  math  model  can  contain  equality 
constraints,  discontinuities  and  break-lines  and  single  point  constraints 
(minimum,  maximum,  saddle  etc.  points).  Furthermore,  the  math  model 
allows  automatic  bridging  of  "smooth  areas"  (sparse  data  sampling)  or 
a  priori  identified  "blunder  areas"  (sampled  data  with  zero  a  priori 
weights).  Those  features  allow  introduction  of  batches  of  fill-in  samples 
replacing  large  areas  of  blunderous  observations,  batches  of  overlapping 
data  samples,  etc.  Thus  the  math  model  can  be  used  for  modeling  even 
"pathologically"  difficult  and  ill-behaving  empirical  functions  with 
proper  computational  efficiency  both  in  the  stages  of  forming  the  data 
base  and  in  the  retrieval  and  usage  of  the  stored  data  base. 

The  above  and  many  other  applications  of  the  sequential  array 
algebra  warrant  detailed  investigations.  For  example  some  carefully 
designed  net  adjustment  problems  of  large  dimensions  are  within  the 
capabilities  of  the  above  array  solution. 
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